home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_11 / ql.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.3 KB  |  74 lines  |  [MATF/MATL]
  1. function D = ql(A,epsilon,show)
  2. % A = ql(A,epsilon,show)
  3. % To find the eigenvalues of a symmetric tridigonal matrix.
  4. % The QL method is employed.
  5. % A is an n x n symmetric tridigonal matrix, input.
  6. % epsilon is the tolerance, input.
  7. % D is the vector of eigenvalues, output.
  8. if nargin==1, show = 0; end
  9. if show==1,
  10.   clc,disp('Currently  A = '),disp(A)
  11.   pause(0.8);
  12. end
  13. [n,n] = size(A);
  14. shift = 0;
  15. m = 1;
  16. cnt = 0;
  17. for m=1:(n-2)  % Loop for eliminating super diagonal.
  18.   cond = 0;
  19.      while cond==0                            % Start of inner iteration loop.
  20.     rts =  eig(A(m:m+1,m:m+1));            % Start FindShift
  21.     r1 = rts(1);                           % Use  eig  to find the
  22.     r2 = rts(2);                           % eigenvalue of A(m:m+1,m:m+1) 
  23.     sh = r1;                               % is closest to A(m,m).
  24.     if (abs(A(m,m)-r2)<abs(A(m,m)-r1)), sh = rts(2); end
  25.     cnt = cnt+1;
  26.     home; if cnt==1, clc; end;
  27.     disp(['The shift is  ',num2str(sh)])   % End of FindShift
  28.        if (abs(A(m,m+1))>epsilon),            % Start of form eigenvalue.
  29.       shift = shift + sh;  
  30.       if show==1,
  31.         disp(''),disp(''),disp('');
  32.       end
  33.       A(m:n,m:n) = A(m:n,m:n) - sh*diag(ones(1,n-m+1));
  34.        else
  35.          A(m,m) = A(m,m) + shift;
  36.       if show==1,
  37.            disp('The eigenvalue is '),disp(A(m,m));
  38.       end
  39.          m = m+1;
  40.          cond = 1;
  41.     end                                    % End of form eigenvalue.
  42.     % Now use the Givens rotations to zero the elements.
  43.     Q = eye(n);
  44.       for j = n-1:-1:m,
  45.         Xp = A(j,j+1);
  46.         Xq = A(j+1,j+1);
  47.         Yq = sqrt(A(j,j+1)^2 +A(j+1,j+1)^2);
  48.         c = -Xq*Yq/(Xp^2+Xq^2);
  49.         s = Xp*Yq/(Xp^2+Xq^2);
  50.         R = [c s; -s c];
  51.         A(j:j+1,:) = R*A(j:j+1,:);
  52.         A(:,j:j+1) = A(:,j:j+1)*R';
  53.       end                                  % End of inner iteration loop
  54.     if show==1,
  55.       Mx = 'QL iteration No. ';
  56.       disp(''),disp([Mx,int2str(cnt)]),disp(''),...
  57.          disp('Currently  A = '),disp(A);
  58.       pause(0.8);
  59.     end
  60.      end
  61. end
  62. % Use  eig  to help find the last two eigenvalues.
  63. A(m:m+1,m:m+1) = shift + diag(eig(A(m:m+1,m:m+1)));
  64. A(n-1,n) = 0;
  65. A(n,n-1) = 0;
  66. if show==1,
  67.   home;disp('');disp('');disp('');disp('');
  68.   Mx = 'QL iteration No. ';
  69.   disp(''),disp([Mx,int2str(cnt)]),disp(''),...
  70.      disp('Currently  A = '),disp(A);
  71.   pause(0.8);
  72. end
  73. D = diag(A);
  74.